The goal here is to use conventional alpha diversity metrics to see how Chao1 richness, shannon diversity and evenness change across samples and to compare those to the values seen using breakaway in the AlphaDiversity.Rmd file
Run AlphaDiversity in scratchnotebooks That file calculates richness in breakawy which I will combine here
#source(here::here("RScripts", "InitialProcessing_3.R"))
source(here::here("RLibraries", "ChesapeakePersonalLibrary.R"))
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
── Attaching packages ──────────────────────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.4.0 ✔ purrr 0.3.4
✔ tibble 3.1.8 ✔ dplyr 1.0.10
✔ tidyr 1.2.1 ✔ stringr 1.4.1
✔ readr 2.1.3 ✔ forcats 0.5.2 ── Conflicts ─────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ksource(here::here("ActiveNotebooks", "BreakawayAlphaDiversity.Rmd"))
processing file: /home/jacob/Projects/ChesapeakeMainstemAnalysis_ToShare/ActiveNotebooks/BreakawayAlphaDiversity.Rmd
|
| | 0%
|
|.... | 4%
|
|......... | 9%
|
|............. | 13%
|
|................. | 17%
|
|...................... | 22%
|
|.......................... | 26%
|
|.............................. | 30%
|
|.................................. | 35%
|
|....................................... | 39%
|
|........................................... | 43%
|
|............................................... | 48%
|
|.................................................... | 52%
|
|........................................................ | 57%
|
|............................................................ | 61%
|
|................................................................. | 65%
|
|..................................................................... | 70%
|
|......................................................................... | 74%
|
|............................................................................. | 78%
|
|.................................................................................. | 83%
|
|...................................................................................... | 87%
|
|.......................................................................................... | 91%
|
|............................................................................................... | 96%
|
|...................................................................................................| 100%
output file: /tmp/RtmpjHvOEV/file136755260521
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 methods overwritten by 'htmltools':
method from
print.html tools:rstudio
print.shiny.tag tools:rstudio
print.shiny.tag.list tools:rstudio
Attaching package: ‘flextable’
The following object is masked from ‘package:purrr’:
compose
Attaching package: ‘ftExtra’
The following object is masked from ‘package:flextable’:
separate_header
Warning: Assuming taxa are rows
library(vegan)
Loading required package: permute
Loading required package: lattice
This is vegan 2.6-3
library(cowplot)
library(flextable)
library(ftExtra)
This file is dedicated to conventional, non div-net/breakaway stats, since breakaway seems to choke on this data.
Reshape back into an ASV matrix, but this time correcting for total abundance
raDf <- nonSpikes_Remake %>% pivot_wider(id_cols = ID, names_from = ASV, values_from = RA, values_fill = 0)
raMat <- raDf %>% column_to_rownames("ID")
raMat1 <- raMat %>% as.matrix()
countMat <- nonSpikes_Remake %>%
pivot_wider(id_cols = ID, names_from = ASV, values_from = reads, values_fill = 0) %>%
column_to_rownames("ID") %>% as.matrix()
seqDep <- countMat %>% apply(1, sum)
min(seqDep)
[1] 852
sampleRichness <- rarefy(countMat, min(seqDep))
rarefy everything to the minimum depth (852)
countRare <- rrarefy(countMat, min(seqDep))
Gamma diversity
specpool(countRare)
Doesn’t finish
#specpool(countMat)
All richness estimates
richnessRare <- estimateR(countRare)
Shannon diversity
shan <- diversity(countRare)
shan
3-1-B-0-2 3-1-B-1-2 3-1-B-180 3-1-B-20 3-1-B-5 3-1-B-500 3-1-B-53 3-1-S-0-2 3-1-S-1-2
4.310891 5.165546 4.689957 5.913745 5.111418 3.870738 5.557505 4.510650 4.848269
3-1-S-180 3-1-S-20 3-1-S-5 3-2-B-0-2 3-2-B-1-2 3-2-B-180 3-2-B-20 3-2-B-5 3-2-B-500
4.696624 5.233180 4.817138 4.421646 4.789616 4.646308 5.183631 5.534260 4.990589
3-2-B-53 3-2-S-0-2 3-2-S-1-2 3-2-S-180 3-2-S-20 3-2-S-5 3-2-S-500 3-2-S-53 3-3-B-0-2
4.976637 3.723508 4.915513 4.789992 4.917679 5.140993 4.974091 4.328113 4.365590
3-3-B-1-2 3-3-B-180 3-3-B-20 3-3-B-5 3-3-B-500 3-3-B-53 3-3-S-180 3-3-S-20 3-3-S-500
4.883612 3.466347 5.650654 5.357062 5.314725 5.519214 4.983904 4.954011 4.912663
3-3-S-53 4-3-B-0-2 4-3-B-1-2 4-3-B-180 4-3-B-20 4-3-B-5 4-3-B-500 4-3-B-53 4-3-O-1-2
4.427651 4.278936 4.953905 4.517047 4.315042 4.678940 4.416019 4.274373 5.007986
4-3-O-180 4-3-O-5 4-3-O-500 4-3-O-53 4-3-S-0-2 4-3-S-180 4-3-S-20 4-3-S-500 4-3-S-53
4.717124 5.223329 4.046068 4.686227 2.805419 4.501172 4.556331 4.744078 4.495019
5-1-S-1-2 5-1-S-180 5-1-S-20 5-1-S-5 5-1-S-500 5-1-S-53 5-5-B-0-2 5-5-B-180 5-5-B-5
4.457156 4.528074 4.165735 3.886173 4.045319 4.008392 4.584748 5.178673 5.438522
5-5-B-500 5-5-B-53 5-5-S-180 5-5-S-5 5-5-S-500 5-5-S-53 C_5P1B_0P2 C_5P1B_180 C_5P1B_1P2
5.030883 4.915445 4.358304 4.955525 4.838466 4.253475 4.239569 4.724150 4.847382
C_5P1B_20 C_5P1B_500 C_5P1B_53
5.376668 4.837291 5.029040
Evenness
pielouJ <- shan/richnessRare["S.chao1",]
pielouJ
3-1-B-0-2 3-1-B-1-2 3-1-B-180 3-1-B-20 3-1-B-5 3-1-B-500 3-1-B-53 3-1-S-0-2 3-1-S-1-2
0.012002661 0.005695199 0.012624379 0.003374655 0.007066477 0.041845812 0.008999821 0.008826596 0.008569315
3-1-S-180 3-1-S-20 3-1-S-5 3-2-B-0-2 3-2-B-1-2 3-2-B-180 3-2-B-20 3-2-B-5 3-2-B-500
0.009531201 0.008122208 0.008588573 0.011376936 0.008519080 0.012306848 0.006157532 0.007217634 0.010603357
3-2-B-53 3-2-S-0-2 3-2-S-1-2 3-2-S-180 3-2-S-20 3-2-S-5 3-2-S-500 3-2-S-53 3-3-B-0-2
0.010199187 0.019991991 0.008528156 0.008806712 0.008659411 0.008230593 0.008635575 0.013009545 0.014582521
3-3-B-1-2 3-3-B-180 3-3-B-20 3-3-B-5 3-3-B-500 3-3-B-53 3-3-S-180 3-3-S-20 3-3-S-500
0.007702219 0.066660513 0.005222085 0.006729824 0.006906726 0.004171370 0.009138070 0.008108037 0.010812581
3-3-S-53 4-3-B-0-2 4-3-B-1-2 4-3-B-180 4-3-B-20 4-3-B-5 4-3-B-500 4-3-B-53 4-3-O-1-2
0.009967698 0.011272718 0.009273930 0.008735626 0.009611597 0.008821531 0.007439386 0.009714484 0.008645639
4-3-O-180 4-3-O-5 4-3-O-500 4-3-O-53 4-3-S-0-2 4-3-S-180 4-3-S-20 4-3-S-500 4-3-S-53
0.010714968 0.008686394 0.015126253 0.010410610 0.124685269 0.014872533 0.010144341 0.009760949 0.012606223
5-1-S-1-2 5-1-S-180 5-1-S-20 5-1-S-5 5-1-S-500 5-1-S-53 5-5-B-0-2 5-5-B-180 5-5-B-5
0.011975164 0.012247374 0.013105682 0.018611254 0.016107958 0.014028890 0.008750354 0.008721128 0.006220032
5-5-B-500 5-5-B-53 5-5-S-180 5-5-S-5 5-5-S-500 5-5-S-53 C_5P1B_0P2 C_5P1B_180 C_5P1B_1P2
0.008597596 0.009145546 0.013668888 0.010845538 0.009159701 0.014943540 0.008794075 0.007531669 0.005760406
C_5P1B_20 C_5P1B_500 C_5P1B_53
0.004266727 0.007981538 0.005106721
diversityData <- sampleData %>%
left_join(richnessRare %>% t() %>% as.data.frame() %>% rownames_to_column("ID"), by = "ID") %>%
left_join(shan %>% enframe(name = "ID", value = "shannonH"), by = "ID") %>%
left_join(pielouJ %>% enframe(name = "ID", value = "pielouJ"), by = "ID") %>%
arrange(Size_Class)
Parameters for all plots
plotSpecs <- list(
facet_wrap(~Depth, ncol = 1) ,
theme_bw(base_size = 16) ,
geom_point(size = 4) ,
geom_path(aes(color = as.factor(Station))) ,
scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
#scale_y_log10nice() ,
scale_shape_manual(values = rep(21:25, 2)) ,
scale_fill_viridis_d(option = "plasma") ,
scale_color_viridis_d(option = "plasma") ,
labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
theme(legend.position = "none",
plot.margin = unit(c(0, 0, 0, 0), "cm"),
axis.title.x = element_blank(),
axis.text.x = element_text(angle = 90, vjust = .5),
axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)
Observed species counts, on rarefied data
plotObs <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = S.obs, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +ylab("Observed ASVs (Rarefied)")#+ scale_y_log10()
plotObs
Estemated Chao1 Richness
plotChao1 <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = S.chao1, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
geom_errorbar(aes(ymin = S.chao1 -2 * se.chao1, ymax = S.chao1 + 2* se.chao1), width = -.1) +
scale_y_log10() +
ylab("Richness (Chao1)")
plotChao1
Shannon diversity
plotShan <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = shannonH, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
ylab("Diversity (Shannon H)") +
lims(y = c(2.5, 6))
plotShan
Evenness
plotPielou <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +scale_y_log10() +ylab("Evenness (PielouJ)")
plotPielou
All plots together
plotAlpha <- plot_grid(plotObs, plotChao1, plotShan, plotPielou, nrow = 1, labels = LETTERS)
plotAlpha
ggsave(here::here("Figures", "ConventionalAlpha.png"), plotAlpha, width = 11, height = 4)
Rarefied observed species numbers
obsMod <- lm(S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(obsMod)
Call:
lm(formula = S.obs ~ log(Size_Class) + I(log(Size_Class)^2) +
I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-222.409 -42.169 3.618 50.716 200.239
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 228869.125 135148.104 1.693 0.094877 .
log(Size_Class) 32.876 8.188 4.015 0.000149 ***
I(log(Size_Class)^2) -6.386 1.523 -4.192 8.07e-05 ***
lat -11916.986 7042.393 -1.692 0.095123 .
I(lat^2) 155.218 91.678 1.693 0.094952 .
depth 3.348 3.243 1.032 0.305574
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 76.98 on 69 degrees of freedom
Multiple R-squared: 0.2574, Adjusted R-squared: 0.2035
F-statistic: 4.782 on 5 and 69 DF, p-value: 0.0008238
Rarified chao1 estimates
chao1Mod <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(chao1Mod)
Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) +
I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-537.58 -134.92 -16.72 110.55 1001.08
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 947118.919 443379.266 2.136 0.036218 *
log(Size_Class) 85.353 26.863 3.177 0.002224 **
I(log(Size_Class)^2) -17.651 4.998 -3.532 0.000741 ***
lat -49369.792 23103.920 -2.137 0.036157 *
I(lat^2) 643.296 300.769 2.139 0.035990 *
depth 19.724 10.640 1.854 0.068034 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 252.5 on 69 degrees of freedom
Multiple R-squared: 0.2128, Adjusted R-squared: 0.1557
F-statistic: 3.729 on 5 and 69 DF, p-value: 0.004785
As above but without latitude and depth
chao1ModSimple <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), data = diversityData)
summary(chao1ModSimple)
Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2),
data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-468.19 -159.00 -14.57 110.38 1095.72
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 561.748 45.961 12.222 < 2e-16 ***
log(Size_Class) 85.808 27.295 3.144 0.002424 **
I(log(Size_Class)^2) -18.065 5.077 -3.558 0.000666 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 256.9 on 72 degrees of freedom
Multiple R-squared: 0.1498, Adjusted R-squared: 0.1262
F-statistic: 6.342 on 2 and 72 DF, p-value: 0.002904
shanMod <- lm(shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(shanMod)
Call:
lm(formula = shannonH ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-1.3141 -0.1915 0.0272 0.3175 0.7540
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.361e+03 7.921e+02 1.718 0.0903 .
log(Size_Class) 2.068e-01 4.799e-02 4.310 5.31e-05 ***
I(log(Size_Class)^2) -3.698e-02 8.929e-03 -4.142 9.60e-05 ***
lat -7.064e+01 4.127e+01 -1.712 0.0915 .
I(lat^2) 9.198e-01 5.373e-01 1.712 0.0914 .
depth 1.177e-02 1.901e-02 0.619 0.5379
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4512 on 69 degrees of freedom
Multiple R-squared: 0.2981, Adjusted R-squared: 0.2472
F-statistic: 5.86 on 5 and 69 DF, p-value: 0.0001437
pielouMod <- lm(pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(pielouMod)
Call:
lm(formula = pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-0.013980 -0.004919 -0.002022 0.000656 0.099590
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.323e+01 2.665e+01 -0.872 0.3865
log(Size_Class) -4.000e-03 1.615e-03 -2.478 0.0157 *
I(log(Size_Class)^2) 6.964e-04 3.004e-04 2.318 0.0234 *
lat 1.209e+00 1.389e+00 0.871 0.3869
I(lat^2) -1.572e-02 1.808e-02 -0.870 0.3875
depth -3.007e-04 6.395e-04 -0.470 0.6397
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.01518 on 69 degrees of freedom
Multiple R-squared: 0.09907, Adjusted R-squared: 0.03379
F-statistic: 1.518 on 5 and 69 DF, p-value: 0.1957
uomisto H (2010a). “A diversity of beta diver- sities: straightening up a concept gone awry. 1. Defining beta diversity as a function of alpha and gamma diversity.” Ecography, 33, 2–2
predict(obsMod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12
226.1310 226.1310 200.6366 200.6366 204.4382 161.4112 161.4112 181.9116 200.7873 301.3648 301.3648 275.8704
13 14 15 16 17 18 19 20 21 22 23 24
275.8704 279.6720 236.6450 236.6450 257.1453 257.1453 331.9534 331.9534 306.4590 306.4590 310.2606 267.2336
25 26 27 28 29 30 31 32 33 34 35 36
267.2336 287.7339 306.6096 306.6096 336.7613 336.7613 311.2669 311.2669 315.0685 315.0685 272.0415 272.0415
37 38 39 40 41 42 43 44 45 46 47 48
292.5418 292.5418 325.4489 299.9545 299.9545 303.7561 303.7561 260.7291 260.7291 260.7291 281.2294 281.2294
49 50 51 52 53 54 55 56 57 58 59 60
300.1051 300.1051 294.1015 294.1015 268.6071 268.6071 272.4087 272.4087 229.3817 229.3817 229.3817 249.8820
61 62 63 64 65 66 67 68 69 70 71 72
249.8820 268.7577 268.7577 253.2659 227.7715 227.7715 231.5731 231.5731 188.5461 188.5461 188.5461 209.0464
73 74 75
209.0464 227.9221 227.9221
$se.fit
[1] 26.93010 26.93010 26.19587 26.19587 26.05607 27.72062 27.72062 29.27714 33.96388 18.89110 18.89110
[12] 18.27607 18.27607 17.16336 19.99663 19.99663 21.45631 21.45631 19.21214 19.21214 18.75132 18.75132
[23] 17.19863 20.04323 20.04323 21.15902 28.11574 28.11574 19.20997 19.20997 18.71830 18.71830 16.92495
[34] 16.92495 19.56830 19.56830 20.53067 20.53067 18.52166 17.88401 17.88401 15.96344 15.96344 18.40305
[45] 18.40305 18.40305 19.39190 19.39190 26.49688 26.49688 19.27023 19.27023 18.38518 18.38518 16.64424
[56] 16.64424 18.38669 18.38669 18.38669 19.45656 19.45656 26.05823 26.05823 24.40667 23.45563 23.45563
[67] 22.29657 22.29657 23.09072 23.09072 23.09072 24.08798 24.08798 29.16339 29.16339
$df
[1] 69
$residual.scale
[1] 76.97868
diversityData$pred_obs = predict(obsMod, se.fit = TRUE)$fit
diversityData$se_obs = predict(obsMod, se.fit = TRUE)$se.fit
plotSpecs2 <- list(
facet_wrap(~Depth, ncol = 1) ,
theme_bw(base_size = 16) ,
#geom_point(size = 4) ,
geom_path(aes(color = as.factor(Station))) ,
scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
#scale_y_log10nice() ,
scale_shape_manual(values = rep(21:25, 2)) ,
scale_fill_viridis_d(option = "plasma") ,
scale_color_viridis_d(option = "plasma") ,
labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
theme(legend.position = "none",
plot.margin = unit(c(0, 0, 0, 0), "cm"),
axis.title.x = element_blank(),
axis.text.x = element_text(angle = 90, vjust = .5),
axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)
plotObs_pred <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_obs, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_obs - 2 * se_obs, yend = pred_obs + 2 * se_obs, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5) +
plotSpecs2 + ylab("Predicted ASVs")
plotObs_pred
predict(chao1Mod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12
470.9464 470.9464 357.1670 357.1670 439.2520 269.8295 269.8295 392.2784 361.0917 669.0136 669.0136 555.2343
13 14 15 16 17 18 19 20 21 22 23 24
555.2343 637.3192 467.8968 467.8968 590.3456 590.3456 745.6874 745.6874 631.9080 631.9080 713.9930 544.5705
25 26 27 28 29 30 31 32 33 34 35 36
544.5705 667.0194 635.8327 635.8327 751.3244 751.3244 637.5450 637.5450 719.6300 719.6300 550.2076 550.2076
37 38 39 40 41 42 43 44 45 46 47 48
672.6564 672.6564 714.6754 600.8960 600.8960 682.9810 682.9810 513.5586 513.5586 513.5586 636.0074 636.0074
49 50 51 52 53 54 55 56 57 58 59 60
604.8207 604.8207 621.2768 621.2768 507.4974 507.4974 589.5824 589.5824 420.1599 420.1599 420.1599 542.6088
61 62 63 64 65 66 67 68 69 70 71 72
542.6088 511.4221 511.4221 502.7609 388.9815 388.9815 471.0664 471.0664 301.6440 301.6440 301.6440 424.0929
73 74 75
424.0929 392.9061 392.9061
$se.fit
[1] 88.34936 88.34936 85.94058 85.94058 85.48195 90.94281 90.94281 96.04928 111.42503 61.97588
[11] 61.97588 59.95815 59.95815 56.30769 65.60278 65.60278 70.39155 70.39155 63.02912 63.02912
[21] 61.51731 61.51731 56.42339 65.75566 65.75566 69.41624 92.23906 92.23906 63.02200 63.02200
[31] 61.40898 61.40898 55.52553 55.52553 64.19755 64.19755 67.35480 67.35480 60.76386 58.67191
[41] 58.67191 52.37114 52.37114 60.37472 60.37472 60.37472 63.61886 63.61886 86.92810 86.92810
[51] 63.21969 63.21969 60.31612 60.31612 54.60461 54.60461 60.32105 60.32105 60.32105 63.83099
[61] 63.83099 85.48902 85.48902 80.07075 76.95068 76.95068 73.14818 73.14818 75.75353 75.75353
[71] 75.75353 79.02524 79.02524 95.67609 95.67609
$df
[1] 69
$residual.scale
[1] 252.5433
diversityData$pred_chao1 = predict(chao1Mod, se.fit = TRUE)$fit
diversityData$se_chao1 = predict(chao1Mod, se.fit = TRUE)$se.fit
plotChao1_pred <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_chao1, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_chao1 - 2 * se_chao1, yend = pred_chao1 + 2 * se_chao1, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5) +
plotSpecs2 + ylab("Predictd Richness (Chao1)") + scale_y_log10()
plotChao1_pred
predict(shanMod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12
4.483447 4.483447 4.343687 4.343687 4.274961 4.008305 4.008305 4.078398 4.341968 4.948626 4.948626 4.808865
13 14 15 16 17 18 19 20 21 22 23 24
4.808865 4.740140 4.473483 4.473483 4.543576 4.543576 5.149242 5.149242 5.009481 5.009481 4.940756 4.674100
25 26 27 28 29 30 31 32 33 34 35 36
4.674100 4.744193 5.007763 5.007763 5.199876 5.199876 5.060115 5.060115 4.991390 4.991390 4.724733 4.724733
37 38 39 40 41 42 43 44 45 46 47 48
4.794826 4.794826 5.150379 5.010618 5.010618 4.941893 4.941893 4.675236 4.675236 4.675236 4.745329 4.745329
49 50 51 52 53 54 55 56 57 58 59 60
5.008900 5.008900 4.988925 4.988925 4.849165 4.849165 4.780439 4.780439 4.513783 4.513783 4.513783 4.583876
61 62 63 64 65 66 67 68 69 70 71 72
4.583876 4.847446 4.847446 4.769215 4.629455 4.629455 4.560730 4.560730 4.294073 4.294073 4.294073 4.364166
73 74 75
4.364166 4.627736 4.627736
$se.fit
[1] 0.15783393 0.15783393 0.15353071 0.15353071 0.15271137 0.16246707 0.16246707 0.17158964 0.19905804
[10] 0.11071836 0.11071836 0.10711375 0.10711375 0.10059229 0.11719773 0.11719773 0.12575275 0.12575275
[19] 0.11259995 0.11259995 0.10989914 0.10989914 0.10079898 0.11747085 0.11747085 0.12401038 0.16478279
[28] 0.16478279 0.11258723 0.11258723 0.10970561 0.10970561 0.09919497 0.09919497 0.11468732 0.11468732
[37] 0.12032767 0.12032767 0.10855312 0.10481590 0.10481590 0.09355973 0.09355973 0.10785794 0.10785794
[46] 0.10785794 0.11365350 0.11365350 0.15529489 0.15529489 0.11294041 0.11294041 0.10775324 0.10775324
[55] 0.09754977 0.09754977 0.10776205 0.10776205 0.10776205 0.11403247 0.11403247 0.15272400 0.15272400
[64] 0.14304440 0.13747047 0.13747047 0.13067740 0.13067740 0.13533179 0.13533179 0.13533179 0.14117663
[73] 0.14117663 0.17092295 0.17092295
$df
[1] 69
$residual.scale
[1] 0.4511623
diversityData$pred_shanH = predict(shanMod, se.fit = TRUE)$fit
diversityData$se_shanH = predict(shanMod, se.fit = TRUE)$se.fit
plotShannonH_pred <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_shanH, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_shanH - 2 * se_shanH, yend = pred_shanH + 2 * se_shanH, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5) +
plotSpecs2 + ylab("Predicted Diversity (Shannon H)") #+ scale_y_log10()
plotShannonH_pred
predict(pielouMod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9
0.019539795 0.019539795 0.021890799 0.021890799 0.021592773 0.025095028 0.025095028 0.022774559 0.018846627
10 11 12 13 14 15 16 17 18
0.010590975 0.010590975 0.012941978 0.012941978 0.012643952 0.016146208 0.016146208 0.013825738 0.013825738
19 20 21 22 23 24 25 26 27
0.006662563 0.006662563 0.009013567 0.009013567 0.008715541 0.012217797 0.012217797 0.009897327 0.005969395
28 29 30 31 32 33 34 35 36
0.005969395 0.005562772 0.005562772 0.007913776 0.007913776 0.007615750 0.007615750 0.011118005 0.011118005
37 38 39 40 41 42 43 44 45
0.008797536 0.008797536 0.006391978 0.008742981 0.008742981 0.008444956 0.008444956 0.011947211 0.011947211
46 47 48 49 50 51 52 53 54
0.011947211 0.009626742 0.009626742 0.005698810 0.005698810 0.009303238 0.009303238 0.011654241 0.011654241
55 56 57 58 59 60 61 62 63
0.011356216 0.011356216 0.014858471 0.014858471 0.014858471 0.012538002 0.012538002 0.008610070 0.008610070
64 65 66 67 68 69 70 71 72
0.013332733 0.015683736 0.015683736 0.015385710 0.015385710 0.018887966 0.018887966 0.018887966 0.016567496
73 74 75
0.016567496 0.012639564 0.012639564
$se.fit
[1] 0.005310081 0.005310081 0.005165306 0.005165306 0.005137741 0.005465956 0.005465956 0.005772871
[9] 0.006697004 0.003724950 0.003724950 0.003603678 0.003603678 0.003384274 0.003942939 0.003942939
[17] 0.004230759 0.004230759 0.003788253 0.003788253 0.003697389 0.003697389 0.003391228 0.003952127
[25] 0.003952127 0.004172140 0.005543865 0.005543865 0.003787825 0.003787825 0.003690877 0.003690877
[33] 0.003337263 0.003337263 0.003858479 0.003858479 0.004048241 0.004048241 0.003652104 0.003526371
[41] 0.003526371 0.003147674 0.003147674 0.003628715 0.003628715 0.003628715 0.003823698 0.003823698
[49] 0.005224659 0.005224659 0.003799707 0.003799707 0.003625193 0.003625193 0.003281913 0.003281913
[57] 0.003625490 0.003625490 0.003625490 0.003836448 0.003836448 0.005138166 0.005138166 0.004812510
[65] 0.004624984 0.004624984 0.004396441 0.004396441 0.004553031 0.004553031 0.004553031 0.004749672
[73] 0.004749672 0.005750441 0.005750441
$df
[1] 69
$residual.scale
[1] 0.01517867
diversityData$pred_pielouJ = predict(pielouMod, se.fit = TRUE)$fit
diversityData$se_pielouJ = predict(pielouMod, se.fit = TRUE)$se.fit
plot_pielouJ_pred <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_pielouJ - 2 * se_pielouJ, yend = pred_pielouJ + 2 * se_pielouJ, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both")) +
plotSpecs2 + ylab("Predicted Evenness (Pielou J)") + scale_y_log10()
plot_pielouJ_pred
plotPredictions <- plot_grid(plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred, nrow = 1, labels = LETTERS)
Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (`geom_segment()`).
plotPredictions
ggsave(here::here("Figures", "ConventionalAlphaPredictions.png"), plotPredictions, width = 11, height = 4)
plot_grid(plotObs, plotChao1, plotShan, plotPielou,
plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred,
nrow = 2, labels = LETTERS)
Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (`geom_segment()`).
alphaSummary <- tibble(
metric = c("Observed ASVs", "Richness (Chao1)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
model = list(obsMod, chao1Mod, shanMod, pielouMod)
)
alphaSummary <- alphaSummary %>%
mutate(df = map(model, ~broom::tidy(summary(.))))
alphaSummary <- alphaSummary %>%
select(-model) %>%
unnest(df)
alphaSummary <- alphaSummary %>%
rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
Term = str_replace(Term, "\\^2", "\\^2\\^"),
Term = str_replace(Term, "depth", "Depth"),
Term = str_replace(Term, "lat", "Latitude"),
Term = str_replace(Term, "_", " ")# BOOKMARK!!
) %>%
mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
rename(`Standard\nError` = `Standard Error`) %>%
identity()
alphaSummary %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>%
bold(i = ~ p< 0.05, j = "p") %>%
colformat_md() %>%
set_table_properties(layout = "autofit") %>%
align(j = -c(1:2), align = "right")
Metric | Term | Estimate | Standard | T Value | p |
Observed ASVs | Intercept | 2.3 x 105 | 1.4 x 105 | 1.69 | 0.095 |
log(Size Class) | 3.3 x 101 | 8.2 x 100 | 4.01 | < 0.001 | |
log(Size Class)2 | -6.4 x 100 | 1.5 x 100 | -4.19 | < 0.001 | |
Latitude | -1.2 x 104 | 7.0 x 103 | -1.69 | 0.095 | |
Latitude2 | 1.6 x 102 | 9.2 x 101 | 1.69 | 0.095 | |
Depth | 3.3 x 100 | 3.2 x 100 | 1.03 | 0.306 | |
Richness (Chao1) | Intercept | 9.5 x 105 | 4.4 x 105 | 2.14 | 0.036 |
log(Size Class) | 8.5 x 101 | 2.7 x 101 | 3.18 | 0.002 | |
log(Size Class)2 | -1.8 x 101 | 5.0 x 100 | -3.53 | < 0.001 | |
Latitude | -4.9 x 104 | 2.3 x 104 | -2.14 | 0.036 | |
Latitude2 | 6.4 x 102 | 3.0 x 102 | 2.14 | 0.036 | |
Depth | 2.0 x 101 | 1.1 x 101 | 1.85 | 0.068 | |
Diversity (Shannon H) | Intercept | 1.4 x 103 | 7.9 x 102 | 1.72 | 0.090 |
log(Size Class) | 2.1 x 10-1 | 4.8 x 10-2 | 4.31 | < 0.001 | |
log(Size Class)2 | -3.7 x 10-2 | 8.9 x 10-3 | -4.14 | < 0.001 | |
Latitude | -7.1 x 101 | 4.1 x 101 | -1.71 | 0.091 | |
Latitude2 | 9.2 x 10-1 | 5.4 x 10-1 | 1.71 | 0.091 | |
Depth | 1.2 x 10-2 | 1.9 x 10-2 | 0.62 | 0.538 | |
Evenness (Pielou J) | Intercept | -2.3 x 101 | 2.7 x 101 | -0.87 | 0.386 |
log(Size Class) | -4.0 x 10-3 | 1.6 x 10-3 | -2.48 | 0.016 | |
log(Size Class)2 | 7.0 x 10-4 | 3.0 x 10-4 | 2.32 | 0.023 | |
Latitude | 1.2 x 100 | 1.4 x 100 | 0.87 | 0.387 | |
Latitude2 | -1.6 x 10-2 | 1.8 x 10-2 | -0.87 | 0.388 | |
Depth | -3.0 x 10-4 | 6.4 x 10-4 | -0.47 | 0.640 |
richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.))))
Warning: `rename_()` was deprecated in dplyr 0.7.0.
Please use `rename()` instead.
diversityDataWB <- full_join(diversityData,
richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.)))),
by = c("ID" = "break_sample_names"), suffix = c("", "_break")) %>%
mutate(pielouJ2 = shannonH/break_estimate) %>%
identity()
diversityDataWB
pielouMod2 <- lm(pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityDataWB)
summary(pielouMod2)
Call:
lm(formula = pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityDataWB)
Residuals:
Min 1Q Median 3Q Max
-0.013954 -0.005126 -0.002510 0.000872 0.105965
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.975e+01 2.757e+01 -0.716 0.4762
log(Size_Class) -3.283e-03 1.670e-03 -1.966 0.0533 .
I(log(Size_Class)^2) 5.746e-04 3.108e-04 1.849 0.0687 .
lat 1.027e+00 1.436e+00 0.715 0.4771
I(lat^2) -1.334e-02 1.870e-02 -0.713 0.4780
depth -2.417e-04 6.615e-04 -0.365 0.7159
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.0157 on 69 degrees of freedom
Multiple R-squared: 0.06698, Adjusted R-squared: -0.0006297
F-statistic: 0.9907 on 5 and 69 DF, p-value: 0.4299
Ok. So the narrative makes sense. Alpha diveristy is driven by variability in richness rather than evenness. Why would we see an effect in chao1 but not breakaway? Because chao1 is more driven by abundant stuff that makes the rarification threshold. My first hunch is that chao1 responds to evenness, but actually that shouldn’t have any effect since there is no evenness variability? Or maybe just that breakaway is more variable (because it detects fine level differences in rare species) and that doesn’t map as nicely with overall patterns.
plotBreak <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
#scale_y_log10()+
ylab("Richness (Breakaway)")
plotBreak
plotPielou2 <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
#scale_y_log10()+
ylab("Evenness (PielouJ)")
plotPielou2
predict(pielouMod2, se.fit = TRUE)
$fit
1 2 3 4 5 6 7
0.0116707844 0.0116707844 0.0135818312 0.0135818312 0.0133535000 0.0160100159 0.0160100159
8 9 10 11 12 13 14
0.0139468521 0.0098844429 0.0043185124 0.0043185124 0.0062295592 0.0062295592 0.0060012280
15 16 17 18 19 20 21
0.0086577439 0.0086577439 0.0065945801 0.0065945801 0.0011020710 0.0011020710 0.0030131178
22 23 24 25 26 27 28
0.0030131178 0.0027847867 0.0054413025 0.0054413025 0.0033781387 -0.0006842705 -0.0006842705
29 30 31 32 33 34 35
0.0002187110 0.0002187110 0.0021297578 0.0021297578 0.0019014267 0.0019014267 0.0045579425
36 37 38 39 40 41 42
0.0045579425 0.0024947787 0.0024947787 0.0009197555 0.0028308023 0.0028308023 0.0026024712
43 44 45 46 47 48 49
0.0026024712 0.0052589870 0.0052589870 0.0052589870 0.0031958232 0.0031958232 -0.0008665860
50 51 52 53 54 55 56
-0.0008665860 0.0033429274 0.0033429274 0.0052539742 0.0052539742 0.0050256431 0.0050256431
57 58 59 60 61 62 63
0.0076821589 0.0076821589 0.0076821589 0.0056189951 0.0056189951 0.0015565859 0.0015565859
64 65 66 67 68 69 70
0.0066852276 0.0085962744 0.0085962744 0.0083679432 0.0083679432 0.0110244591 0.0110244591
71 72 73 74 75
0.0110244591 0.0089612953 0.0089612953 0.0048988861 0.0048988861
$se.fit
[1] 0.005493164 0.005493164 0.005343397 0.005343397 0.005314882 0.005654413 0.005654413 0.005971910
[9] 0.006927905 0.003853380 0.003853380 0.003727927 0.003727927 0.003500958 0.004078885 0.004078885
[17] 0.004376629 0.004376629 0.003918866 0.003918866 0.003824869 0.003824869 0.003508152 0.004088390
[25] 0.004088390 0.004315988 0.005735008 0.005735008 0.003918423 0.003918423 0.003818133 0.003818133
[33] 0.003452327 0.003452327 0.003991513 0.003991513 0.004187817 0.004187817 0.003778022 0.003647954
[41] 0.003647954 0.003256201 0.003256201 0.003753828 0.003753828 0.003753828 0.003955533 0.003955533
[49] 0.005404797 0.005404797 0.003930715 0.003930715 0.003750184 0.003750184 0.003395068 0.003395068
[57] 0.003750491 0.003750491 0.003750491 0.003968722 0.003968722 0.005315321 0.005315321 0.004978438
[65] 0.004784446 0.004784446 0.004548023 0.004548023 0.004710012 0.004710012 0.004710012 0.004913433
[73] 0.004913433 0.005948707 0.005948707
$df
[1] 69
$residual.scale
[1] 0.015702
diversityDataWB$pred_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$fit
diversityDataWB$se_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$se.fit
plot_pielouJ2_pred <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_pielouJ2 - 2 * se_pielouJ2, yend = pred_pielouJ2 + 2 * se_pielouJ2, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both")) +
plotSpecs2 + ylab("Predicted Evenness (Pielou J2)") #+ scale_y_log10()
plot_pielouJ2_pred
plotBreakaway <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
geom_errorbar(aes(ymin = break_lower, ymax = break_upper), width = -.1) +
scale_y_log10() +
ylab("Richness (Breakaway)")
plotBreakaway
#predict(breakLm, se.fit = TRUE)
# doesn't work because built with a different data frame
Why are these not smooth curves?!! What if I redo the model, this time with the same data frame
breakLm2 <- lm(break_estimate ~ log(Size_Class) + I(log(Size_Class) ^2) + lat + I(lat^2) + depth ,data = diversityDataWB)
breakLm2 %>% summary()
Call:
lm(formula = break_estimate ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityDataWB)
Residuals:
Min 1Q Median 3Q Max
-2974.5 -1191.2 -151.6 599.9 6768.1
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7124615.61 3339862.88 2.133 0.0365 *
log(Size_Class) 244.45 202.35 1.208 0.2312
I(log(Size_Class)^2) -75.16 37.65 -1.996 0.0498 *
lat -370568.38 174035.93 -2.129 0.0368 *
I(lat^2) 4817.28 2265.61 2.126 0.0371 *
depth 151.10 80.15 1.885 0.0636 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1902 on 69 degrees of freedom
Multiple R-squared: 0.1414, Adjusted R-squared: 0.0792
F-statistic: 2.273 on 5 and 69 DF, p-value: 0.0567
Note the non statistical significance overall
#predict(breakLm2, se.fit = TRUE)
diversityDataWB$pred_break = predict(breakLm2, se.fit = TRUE)$fit
diversityDataWB$se_break = predict(breakLm2, se.fit = TRUE)$se.fit
plot_break_pred <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
# filter(Station == 4.3) %>%
ggplot(aes(x = Size_Class, y = pred_break, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_break - 2 * se_break, yend = pred_break + 2 * se_break, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both")) +
plotSpecs2 + ylab("Predicted Richness (Breakaway -- LM)") #+ scale_y_log10()
plot_break_pred
plotAlphaWB <- plot_grid(plotBreakaway, plotShan, plotPielou2, nrow = 1, labels = LETTERS)
plotAlphaWB
ggsave(here::here("Figures", "BreakawayAlpha.png"), plotAlpha, width = 11, height = 4)
Summary table I want both breakaway metrics here
bettaTable <- myBet$table %>%
as.data.frame() %>%
rename(estimate = Estimates,
`std.error` = `Standard Errors`,
`p.value`=`p-values`
) %>%
mutate(`statistic` = NA) %>%
rownames_to_column(var = "term") %>%
select(term, estimate, std.error, statistic, p.value) %>%
as_tibble()
bettaTable
alphaSummary2 <- tibble(
metric = c("Richness (Breakaway -- LM)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
model = list(breakLm, shanMod, pielouMod2)
)
alphaSummary2 <- alphaSummary2 %>%
mutate(df = map(model, ~broom::tidy(summary(.))))
## Add in willis variables
breakawaySummary <- tibble(
metric = "Richness (Breakaway -- Betta)",
model = NULL,
df = list(bettaTable)
)
alphaSummary2 = bind_rows(breakawaySummary, alphaSummary2)
alphaSummary2 <- alphaSummary2 %>%
select(-model) %>%
unnest(df)
alphaSummary2 <- alphaSummary2 %>%
rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
Term = str_replace(Term, "\\^2", "\\^2\\^"),
Term = str_replace(Term, "depth", "Depth"),
Term = str_replace(Term, "lat", "Latitude"),
Term = str_replace(Term, "_", " ")# BOOKMARK!!
) %>%
mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
rename(`Standard\nError` = `Standard Error`) %>%
identity()
alphaSummary2
alphaTable2 <- alphaSummary2 %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>% bold(i = ~ p< 0.05, j = "p") %>% colformat_md() %>% set_table_properties(layout = "autofit") %>%
align(j = -c(1:2), align = "right")
alphaTable2
Metric | Term | Estimate | Standard | T Value | p |
Richness (Breakaway – Betta) | Intercept | 7.1 x 106 | 2.4 x 102 | NA | < 0.001 |
log(Size Class) | 1.2 x 102 | 6.1 x 101 | NA | 0.058 | |
log(Size Class)2 | -5.0 x 101 | 1.2 x 101 | NA | < 0.001 | |
Latitude | -3.7 x 105 | 6.1 x 100 | NA | < 0.001 | |
Latitude2 | 4.8 x 103 | 1.6 x 10-1 | NA | < 0.001 | |
Depth | 1.5 x 102 | 1.0 x 101 | NA | < 0.001 | |
Richness (Breakaway – LM) | Intercept | 7.1 x 106 | 3.3 x 106 | 2.13 | 0.036 |
log(Size Class) | 2.4 x 102 | 2.0 x 102 | 1.21 | 0.231 | |
log(Size Class)2 | -7.5 x 101 | 3.8 x 101 | -2.00 | 0.050 | |
Latitude | -3.7 x 105 | 1.7 x 105 | -2.13 | 0.037 | |
Latitude2 | 4.8 x 103 | 2.3 x 103 | 2.13 | 0.037 | |
Depth | 1.5 x 102 | 8.0 x 101 | 1.89 | 0.064 | |
Diversity (Shannon H) | Intercept | 1.4 x 103 | 7.9 x 102 | 1.72 | 0.090 |
log(Size Class) | 2.1 x 10-1 | 4.8 x 10-2 | 4.31 | < 0.001 | |
log(Size Class)2 | -3.7 x 10-2 | 8.9 x 10-3 | -4.14 | < 0.001 | |
Latitude | -7.1 x 101 | 4.1 x 101 | -1.71 | 0.091 | |
Latitude2 | 9.2 x 10-1 | 5.4 x 10-1 | 1.71 | 0.091 | |
Depth | 1.2 x 10-2 | 1.9 x 10-2 | 0.62 | 0.538 | |
Evenness (Pielou J) | Intercept | -2.0 x 101 | 2.8 x 101 | -0.72 | 0.476 |
log(Size Class) | -3.3 x 10-3 | 1.7 x 10-3 | -1.97 | 0.053 | |
log(Size Class)2 | 5.7 x 10-4 | 3.1 x 10-4 | 1.85 | 0.069 | |
Latitude | 1.0 x 100 | 1.4 x 100 | 0.71 | 0.477 | |
Latitude2 | -1.3 x 10-2 | 1.9 x 10-2 | -0.71 | 0.478 | |
Depth | -2.4 x 10-4 | 6.6 x 10-4 | -0.37 | 0.716 |
alphaTable2 %>% save_as_docx(path = here::here("Tables", "alphaTable2.docx"))
myBet$table
plot_grid(plot_break_pred,plotShannonH_pred,plot_pielouJ2_pred, nrow = 1, labels = LETTERS)